#    pythonequations is a collection of equations expressed as Python classes
#    Copyright (C) 2008 James R. Phillips
#    2548 Vera Cruz Drive
#    Birmingham, AL 35235 USA
#    email: zunzun@zunzun.com
#
#    License: BSD-style (see LICENSE.txt in main source directory)
#    Version info: $Id: __init__.py 267 2010-09-25 13:25:43Z zunzun.com $

import pythonequations, pythonequations.EquationBaseClasses, pythonequations.ExtraCodeForEquationBaseClasses
import numpy, scipy.interpolate
numpy.seterr(all = 'raise') # numpy raises warnings, convert to exceptions to trap them


class Spline2D(pythonequations.EquationBaseClasses.Equation2D):
    RequiresAutoGeneratedGrowthAndDecayForms = False
    RequiresAutoGeneratedOffsetForm = False
    RequiresAutoGeneratedReciprocalForm = False
    RequiresAutoGeneratedInverseForms = False
    _name = "Spline"
    _HTML = "y = B-Spline Interpolation Curve"
    splineFlag = 1
    function_cpp_code = ';' # unused

    def CreateCacheGenerationList(self):
        self.CacheGenerationList = []
        self.CacheGenerationList.append([pythonequations.ExtraCodeForEquationBaseClasses.CG_X(NameOrValueFlag=1), []])


    def EvaluateCachedData(self, coeff, _id):
        return self.scipySpline(_id[0]) # scipySpline should already contains coefficients internally, along with knot points


    def FitToCacheData(self):
        self.FindOrCreateCache()
        
        # sorting along the X axis yields improved results over those from using unsorted data
        data = numpy.array([self.cache[0], self.DependentDataArray])
        data = data.take(numpy.argsort(data)[0], axis=1)
        self.scipySpline = scipy.interpolate.fitpack2.UnivariateSpline(data[0], data[1], s=self.smoothingFactor, k=self.order)
        self.coefficientArray = self.scipySpline.get_coeffs() # scipySpline now has spline coeffs internally, also make them available as "coefficientArray"


    def CalculateFittingTarget(self, in_coeffArray):
        raise NotImplementedError, 'Not implemented for spline interpolations'


    def CodePYTHON(self):
        s  = '''# To the best of my knowledge this code is correct.
# If you find any errors or problems please contact
# me at zunzun@zunzun.com.
#      James
#
# the code below was partially based on the fortran code at:
# http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/splev.f
# http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/fpbspl.f

'''
        s += "def " + self.__class__.__name__ + "_evaluation(x_in):\n"
        s += '    t = ['
        for i in range(len(self.scipySpline._eval_args[0])):
            s += "%-.16E" % (self.scipySpline._eval_args[0][i])
            if i < (len(self.scipySpline._eval_args[0]) - 1):
                s += ", "
        s += ']\n'
        s += '    coeff = ['
        for i in range(len(self.scipySpline._eval_args[1])):
            s += "%-.16E" % (self.scipySpline._eval_args[1][i])
            if i < (len(self.scipySpline._eval_args[1]) - 1):
                s += ", "
        s += ']\n'
        s += '    n = ' + str(len(self.scipySpline._eval_args[0])) + '\n'
        s += '    k = ' + str(self.order) + '\n'
        s += '''
    h = [0.0] * 25

    k1 = k+1
    l = k1
    l1 = l+1

    while x_in < t[l-1] and l1 != (k1+1):
        l1 = l
        l = l-1

    while x_in >= t[l1-1] and l != (n-k1):
        l = l1
        l1 += 1

    hh = [0.0] * 25

    h[0] = 1.0
    for j in range(1, k+1):
        for i in range(j):
            hh[i] = h[i]
        h[0] = 0.0
        for i in range(j):
            li = l+i
            lj = li-j
            if t[li] != t[lj]:
                f = hh[i] / (t[li] - t[lj])
                h[i] = h[i] + f * (t[li] - x_in)
                h[i+1] = f * (x_in - t[lj])
            else:
                h[i+1] = 0.0

    temp = 0.0
    ll = l - k1
    for j in range(k1):
        ll = ll + 1
        temp = temp + coeff[ll-1] * h[j]
        
    return temp
'''
        return s


    def CodeCPP(self):
        s  = '''// To the best of my knowledge this code is correct.
// If you find any errors or problems please contact
// me at zunzun@zunzun.com.
//      James
//
// the code below was partially based on the fortran code at:
// http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/splev.f
// http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/fpbspl.f

'''
        s += "double " + self.__class__.__name__ + "_evaluation(double x_in)\n"
        s += '{\n'
        s += '    double t [] = {'
        for i in range(len(self.scipySpline._eval_args[0])):
            s += "%-.16E" % (self.scipySpline._eval_args[0][i])
            if i < (len(self.scipySpline._eval_args[0]) - 1):
                s += ", "
        s += '};\n'
        s += '    double coeff [] = {'
        for i in range(len(self.scipySpline._eval_args[1])):
            s += "%-.16E" % (self.scipySpline._eval_args[1][i])
            if i < (len(self.scipySpline._eval_args[1]) - 1):
                s += ", "
        s += '};\n'
        s += '    int n = ' + str(len(self.scipySpline._eval_args[0])) + ';\n'
        s += '    int k = ' + str(self.order) + ';\n'
        s += '''
    double h [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    double hh [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

    int i, j, li, lj, ll;
    double f, temp;
    
    int k1 = k+1;
    int l = k1;
    int l1 = l+1;

    while ((x_in < t[l-1]) && (l1 != (k1+1)))
    {
        l1 = l;
        l = l-1;
    }

    while ((x_in >= t[l1-1]) && (l != (n-k1)))
    {
        l = l1;
        l1 += 1;
    }

    h[0] = 1.0;
    for (j = 1; j < k+1; j++)
    {
        for (i = 0; i < j; i++)
        {
            hh[i] = h[i];
        }
        h[0] = 0.0;
        for (i = 0; i < j; i++)
        {
            li = l+i;
            lj = li-j;
            if (t[li] != t[lj])
            {
                f = hh[i] / (t[li] - t[lj]);
                h[i] = h[i] + f * (t[li] - x_in);
                h[i+1] = f * (x_in - t[lj]);
            }
            else
            {
                h[i+1] = 0.0;
            }
        }
    }

    temp = 0.0;
    ll = l - k1;
    for (j = 0; j < k1; j++)
    {
        ll = ll + 1;
        temp = temp + coeff[ll-1] * h[j];
    }
    
    return temp;
}
'''
        return s


    def CodeJAVA(self):
        s  = '''// To the best of my knowledge this code is correct.
// If you find any errors or problems please contact
// me at zunzun@zunzun.com.
//      James
//
// the code below was partially based on the fortran code at:
// http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/splev.f
// http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/fpbspl.f

'''
        s += "class " + self.__class__.__name__ + "\n"
        s += '{\n'
        s += "    double " + self.__class__.__name__ + "_evaluation(double x_in)\n"
        s += '    {\n'
        s += '        double t [] = {'
        for i in range(len(self.scipySpline._eval_args[0])):
            s += "%-.16E" % (self.scipySpline._eval_args[0][i])
            if i < (len(self.scipySpline._eval_args[0]) - 1):
                s += ", "
        s += '};\n'
        s += '        double coeff [] = {'
        for i in range(len(self.scipySpline._eval_args[1])):
            s += "%-.16E" % (self.scipySpline._eval_args[1][i])
            if i < (len(self.scipySpline._eval_args[1]) - 1):
                s += ", "
        s += '};\n'
        s += '        int n = ' + str(len(self.scipySpline._eval_args[0])) + ';\n'
        s += '        int k = ' + str(self.order) + ';\n'
        s += '''
        double h [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
        double hh [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

        int i, j, li, lj, ll;
        double f, temp;
        
        int k1 = k+1;
        int l = k1;
        int l1 = l+1;

        while ((x_in < t[l-1]) && (l1 != (k1+1)))
        {
            l1 = l;
            l = l-1;
        }

        while ((x_in >= t[l1-1]) && (l != (n-k1)))
        {
            l = l1;
            l1 += 1;
        }

        h[0] = 1.0;
        for (j = 1; j < k+1; j++)
        {
            for (i = 0; i < j; i++)
            {
                hh[i] = h[i];
            }
            h[0] = 0.0;
            for (i = 0; i < j; i++)
            {
                li = l+i;
                lj = li-j;
                if (t[li] != t[lj])
                {
                    f = hh[i] / (t[li] - t[lj]);
                    h[i] = h[i] + f * (t[li] - x_in);
                    h[i+1] = f * (x_in - t[lj]);
                }
                else
                {
                    h[i+1] = 0.0;
                }
            }
        }

        temp = 0.0;
        ll = l - k1;
        for (j = 0; j < k1; j++)
        {
            ll = ll + 1;
            temp = temp + coeff[ll-1] * h[j];
        }
        
        return temp;
    }
}
'''
        return s


    def CodeCS(self):
        raise NotImplementedError, 'Not implemented for spline interpolations'


    def CodeSCILAB(self):
        raise NotImplementedError, 'Not implemented for spline interpolations'


    def CodeMATLAB(self):
        raise NotImplementedError, 'Not implemented for spline interpolations'

    
    def CodeVBA(self):
        raise NotImplementedError, 'Not implemented for spline interpolations'
